ERHS732 final project
Grace Kuiper
12/08/2021
Background
Data are publicly available through Dryad for “Stronger social bonds do not always predict greater longevity in a gregarious primate” publication (Thompson and Cords, 2018). For this publication, Thompson and Cords conducted survival analyses to understand the association between measures of sociability and mortality among blue monkeys in Kenya. Using these data, the analysis described here will attempt to answer two follow-up research questions which expand on the original results presented by Thompson and Cords.
Introduction
In this investigation, blue monkey groups were observed for 9 years to understand the impact of sociability on survival. Among humans and non-human species that demonstrate strong social ties, bonding can be an important factor for longevity, offspring survival, and reproductive rate. Together, these are indicators of fitness.
Social measures collected by Thompson and Cords that are available in their publicly available dataset include bond strength and partner consistency. These were evaluated among each female dyad on an annual basis. For each subject, the bond strength and partner consistency scores for her top three and top six closest partners were averaged. Other available demographic and environmental factors include each subject’s start age and age at censorship, mortality outcome, age at first reproduction, dominance rank, and the number of female groupmates.
For time-varying data, subjects’ social and environmental factors were annually for each study year during which they were observed:
| Subject | Age at year start | Age at year end or censorship | Death | Bond strength - 3 closest partners | Bond strength - 6 closest partners | Dominance rank | Number of groupmates | Consistency - 3 closest partners | Consistency - 6 closest partners |
|---|---|---|---|---|---|---|---|---|---|
| adel | 4644 | 5009 | 0 | 3.917304 | 2.513434 | 0.8666667 | 16 | 0.3333333 | 0.5000000 |
| adel | 5010 | 5374 | 0 | 3.603165 | 2.599880 | 0.6250000 | 9 | 0.6666667 | 0.8333333 |
| adel | 5383 | 5739 | 0 | 6.623989 | 3.686815 | 0.6666667 | 10 | 0.3333333 | 0.5000000 |
| adel | 5740 | 5801 | 1 | 2.207133 | 1.241564 | 0.8181818 | 10 | 0.6666667 | 0.6666667 |
| agne | 3673 | 4038 | 0 | 9.184368 | 5.301906 | 0.6000000 | 16 | 0.6666667 | 0.5000000 |
| agne | 4039 | 4403 | 0 | 2.217799 | 1.353885 | 0.2500000 | 9 | 0.3333333 | 0.6666667 |
In their analysis, Thompson and Cords categorized individuals into one of four groups based on both their bond strength and partner consistency: +/+ (bond strength and consistency above the mean), +/- (bond strength above the mean, consistency below the mean), -/+, and -/-. In both same-year and multi-year models, belonging to the +/- group was found to be the most hazardous for mortality (multi-year HR [referent = +/+], 95% CI: 20.3, 3.4-123). The authors hypothesize that these findings indicate a costliness for female blue monkeys to be part of a strongly-bonded partnership. Females that become highly invested in an inconsistent dyad are at a disadvantage.
Research question
Does a decrease in partnership consistency over the previous year affect the likelihood of mortality for female blue monkeys?
In light of Thompson and Cords’ findings, do individuals who experience a change in closest partnerships, as indicated by a decrease in consistency score from lagges years, experience increased mortality risk as compared to those whose partnership consistency stays the same or increases.
Methods
Simple approach
Here the time-varying dataset will be used. A new variable will be created which indicates a decrease in partnership consistency from the previous year versus partnership consistency that stayed the same or increased. This new indicator variable will be added to the model fit by Thompson and Cords, which included each subject’s social category (+/+, +/-, -/+, or -/-), the number of female groupmates, and their dominance rank during each year of observation.
coxph_mod1 <- coxph(Surv(start.age, stop.age, death) ~ as.factor(soc.cat) + no.groupmates + rank +
as.factor(cons.decrease),
data=annual_df)Time-varying IPW approach
The association between potential confounders and the exposure of interest (loss in partner consistency) will be evaluated. Those that are imbalanced across exposure groups, as well as age at study entry as a time-invariant covariate will be used to predict propensity scores of rloss in partner consistency. Cumulative product inverse probability weights will be calculated for each subject. If necessary, a numerator model will be fit to stabilize the IPWs. A Cox proportional hazards model will be fit with any covariates included in the numerator model, as well as the consistency loss exposure variable. This model will include the cumulative product IPW as the weights option and subject ID for a cluster effect.
Results
Data exploration
## There were 20 deaths among the cohort over the course of this study.
Simple approach
A new variable was be created to indicate a decrease in partnership consistency from the previous year versus partnership consistency that stayed the same or increased. This new indicator variable was be added to the model fit by Thompson and Cords.
coxph_mod1 <- coxph(Surv(time, time2, death) ~
a.st.co3 + af.groupmates + rank +
entry.years + cons.p3_loss,
data=annual_df_ps)## # A tibble: 7 x 4
## term hr low_hr high_hr
## <chr> <dbl> <dbl> <dbl>
## 1 a.st.co31 1.25 0.209 7.50
## 2 a.st.co32 2.18 0.515 9.27
## 3 a.st.co33 1.96 0.320 12.0
## 4 af.groupmates 0.965 0.864 1.08
## 5 rank 1.27 0.263 6.12
## 6 entry.years 1.19 1.06 1.34
## 7 cons.p3_loss1 1.54 0.421 5.64## chisq df p
## a.st.co3 2.62183 3 0.45
## af.groupmates 0.00601 1 0.94
## rank 0.01533 1 0.90
## entry.years 1.22770 1 0.27
## cons.p3_loss 0.11903 1 0.73
## GLOBAL 4.03297 7 0.78Time-varying IPW approach
Let’s start by checking whether number of groupmates and dominance ranke are, in fact, associated with exposures (partnership consistency and bond strength).
Are the confounders associated with the exposures?
## # A tibble: 4 x 2
## name value
## <chr> <dbl>
## 1 mean_rank_loss 0.573
## 2 mean_rank_noloss 0.487
## 3 sd_rank 0.323
## 4 stand_diff_rank 0.268Since the standardized mean difference is >0.1, there is an imbalance in rank across the exposure groups (loss vs. no loss in partner consistency).
## # A tibble: 4 x 2
## name value
## <chr> <dbl>
## 1 mean_groupmates_loss 14.2
## 2 mean_groupmates_noloss 13.8
## 3 sd_groupmates 4.79
## 4 stand_diff_groupmates 0.0799So number of groupmates is not imbalanced across exposure groups.
## # A tibble: 2 x 9
## cons.p3_loss n_negneg n_posneg n_negpos n_pospos perc_negneg perc_posneg
## <fct> <int> <int> <int> <int> <dbl> <dbl>
## 1 0 36 70 16 64 0.194 0.376
## 2 1 49 9 24 3 0.576 0.106
## # ... with 2 more variables: perc_negpos <dbl>, perc_pospos <dbl>The social categories are very imbalanced across exposure groups.
Estimating propensity scores
model_ps <- glm(cons.p3_loss ~ ns(time,df=2) + entry.years + rank_l1 + af.groupmates_l1 +
a.st.co3_l1,
family = "binomial", data = annual_df_ps)## # A tibble: 9 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -1.42 0.925 -1.53 0.126
## 2 ns(time, df = 2)1 2.18 1.08 2.03 0.0428
## 3 ns(time, df = 2)2 -0.502 0.556 -0.904 0.366
## 4 entry.years -0.0648 0.0459 -1.41 0.158
## 5 rank_l1 0.331 0.514 0.643 0.520
## 6 af.groupmates_l1 0.0764 0.0387 1.98 0.0482
## 7 a.st.co3_l11 -2.77 0.511 -5.42 0.0000000611
## 8 a.st.co3_l12 0.759 0.383 1.98 0.0474
## 9 a.st.co3_l13 -2.73 0.653 -4.18 0.0000288Calculating IPWs
## # A tibble: 1 x 4
## `Mean cumulative p~ `Minimum cumulativ~ `Maximum cumulativ~ `Sum cumulative p~
## <dbl> <dbl> <dbl> <dbl>
## 1 6.44 1.02 74.0 1744.Fitting numerator model and calculate stabilized IPW
mod_IPWnum_2 <- glm(cons.p3_loss ~ ns(time, df = 2) + entry.years +a.st.co3_l1,
family = "binomial", data = annual_df_ps)## # A tibble: 7 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -0.114 0.653 -0.175 0.861
## 2 ns(time, df = 2)1 2.00 1.06 1.88 0.0595
## 3 ns(time, df = 2)2 -0.382 0.547 -0.698 0.485
## 4 entry.years -0.0727 0.0452 -1.61 0.107
## 5 a.st.co3_l11 -2.57 0.493 -5.22 0.000000183
## 6 a.st.co3_l12 0.825 0.377 2.19 0.0285
## 7 a.st.co3_l13 -2.71 0.648 -4.17 0.0000299## # A tibble: 1 x 4
## `Mean stabilized IPW` `Minimum stabilize~ `Maximum stabiliz~ `Sum stabilized ~
## <dbl> <dbl> <dbl> <dbl>
## 1 0.992 0.502 1.71 269.Cox proportional hazards model
coxph_modIPW_tv <- coxph(Surv(time,time2,death) ~
cons.p3_loss + entry.years + a.st.co3_l1,
weights=w_new_2, cluster=subj,
data=annual_df_ps)## # A tibble: 5 x 6
## term estimate std.error robust.se statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 cons.p3_loss1 0.629 0.651 0.642 0.981 0.327
## 2 entry.years 0.191 0.0545 0.0477 4.00 0.0000628
## 3 a.st.co3_l11 0.787 0.757 0.654 1.20 0.229
## 4 a.st.co3_l12 0.600 0.726 0.691 0.868 0.385
## 5 a.st.co3_l13 0.779 0.843 0.721 1.08 0.280## # A tibble: 5 x 4
## term hr low_hr high_hr
## <chr> <dbl> <dbl> <dbl>
## 1 cons.p3_loss1 1.88 0.534 6.60
## 2 entry.years 1.21 1.10 1.33
## 3 a.st.co3_l11 2.20 0.609 7.92
## 4 a.st.co3_l12 1.82 0.470 7.06
## 5 a.st.co3_l13 2.18 0.531 8.95## chisq df p
## cons.p3_loss 0.488 1 0.48
## entry.years 1.254 1 0.26
## a.st.co3_l1 1.111 3 0.77
## GLOBAL 4.006 5 0.55“The correct interpretation of this HR is as the effect of exposure in the exposed. In other words it is the HR comparing what would have happened if the exposed were exposed (what actually happened) to what would have happened if the exposed were unexposed. By contrast remember that the interpretation of the HR in the case of IPW was the HR comparing what would have happened if everyone was exposed to what would have happened if nobody was exposed.”